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Abstract 

In this paper we propose the first better than second order accurate method in space 
and time for the numerical solution of the resistive relativistic magnetohydrodynam- 
ics (RRMHD) equations on unstructured meshes in multiple space dimensions. The 
nonlinear system under consideration is purely hyperbolic and contains a source 
term, the one for the evolution of the electric field, that becomes stiff for low values 
of the resistivity. 

For the spatial discretization we propose to use high order PnPm schemes as intro- 
duced in [10] for hyperbolic conservation laws and a high order accurate unsplit time 
discretization is achieved using the element-local space-time discontinuous Galerkin 
approach proposed in [11] for one-dimensional balance laws with stiff source terms. 
The divergence free character of the magnetic field is accounted for through the 
divergence cleaning procedure of Dedner et al. [7]. 

To validate our high order method we first solve some numerical test cases for 
which exact analytical reference solutions are known and we also show numerical 
convergence studies in the stiff limit of the RRMHD equations using P^Pm schemes 
from third to fifth order of accuracy in space and time. We also present some 
applications with shock waves such as a classical shock tube problem with different 
values for the conductivity as well as a relativistic MHD rotor problem and the 
relativistic equivalent of the Orszag-Tang vortex problem. We have verified that the 
proposed method can handle equally well the resistive regime and the stiff limit of 
ideal relativistic MHD. For these reasons it provides a powerful tool for relativistic 
astrophysical simulations involving the appearance of magnetic reconnection. 

Key words: resistive relativistic magnetohydrodynamics, unstructured meshes, 
stiff source terms, high order finite volume and discontinuous Galerkin methods, 
PnPm schemes 
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1 Introduction 



Although the assumption of infinite conductivity is often justified in astro- 
physics, there are nevertheless situations in which neglecting the resistivity 
of the plasma may lead to rather inaccurate or simply wrong conclusions. 
This is particularly the case for those physical systems involving processes 
that present magnetic reconnection, such as in the magnetospheres of pul- 
sars near the Y-point, where the outermost magnetic field lines intersect the 
equatorial plane [H], [33]; or in soft gamma-ray repeaters where giant flares 
could be the explanation of the observed strongly magnetized and relativis- 
tic ejection events [25]; or in extragalactic jets, where particle acceleration by 
magnetic reconnection in electron-positron plasmas is supposed to take place 
[IE], PS], and also in active galactic nuclei, where Petschek magnetic recon- 
nection, associated with MHD turbulence, may generate violent releases of 
energy [9]. Moreover, the presence of relativistic motion makes resistive ef- 
fects quantitatively and qualitatively different from those encountered in the 
Newtonian regime. For example, the relativistic reconnection of Petschek type 
with non-strictly parallel reconnecting fields produces a strong compression of 
the plasma and the energy of the reconnecting field can be largely propagated 
away in the form of a Poynting flux [21]. In addition, the reconnection rate is 
also affected, and it is roughly obtained by replacing the Alfven wave velocity 
with the speed of light in the corresponding formulas. 

For all these reasons, and moreover because the question on whether relativis- 
tic magnetic reconnection is an efficient energy converter is still under debate 
(see the discussion in [26] and [36]) there is a strong interest in the numerical 
solution of the full system of RRMHD, by providing a single computational 
tool that can equally well handle situations of low and high resistivity, as com- 
monly encountered in all realistic physical scenarios. The equations to solve 
are particularly challenging from the numerical point of view, since, as recently 
shown by [22] and [28], they become stiff for large values of the conductivity. 
Namely, the RRMHD equations can be cast into the following general form of 
a hyperbolic balance law 

^-W + V-F(W) = S(W), (1) 

where W is the state vector, F(W) is a nonlinear flux tensor that depends on 
the state W and S(W) is a nonlinear source term that becomes stiff at high 
conductivities. In this paper we solve the RRMHD equations by applying the 
high order accurate method recently proposed by [TT] to cope with stiff source 
terms on the right hand side of ([T| and maintaining at the same time better 
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than second order of accuracy in space and time. The numerical method is 
formulated as one-step local predictor global corrector method. The predictor 
is based on an element-local weak solution of ([I]), where inside each element 
the governing PDE ([I]) is solved in the small (see [17]) by means of a locally 
implicit space-time discontinuous Galerkin scheme. This leads to an algebraic 
system of non-linear equations that must be solved individually for each el- 
ement. The globally explicit update in time, on the other hand, is obtained 
by either standard finite volume or discontinuous Galerkin methods, or, fi- 
nally, by a recently proposed generalization of the two named PnPm schemes 
according to [TO] . 

The plan of the paper is as follows. In Sect. [2]we briefly review the peculiar fea- 
tures of the RRMHD equations. The core of the numerical method is described 
in Sect. [3] while Sect. [4] is devoted to the validation of the scheme through a 
large class of numerical tests. Finally, the conclusions are reported in Sect. [5] 
We have considered only flat spacetimes in Cartesian coordinates, namely the 
metric rj^ = diag(— 1, 1, 1, 1), where from now onwards we agree to use Greek 
letters fi, u, A, . . . (running from to 3) for indices of four-dimensional space- 
time tensors, while using Latin letters i,j,k,... (running from 1 to 3) for 
indices of three-dimensional spatial tensors. Finally, we set the speed of light 
c = 1 and make use of the Lorentz-Heaviside notation for the electromagnetic 
quantities, such that all \fkn factors disappear. We use Einstein summation 
convention over repeated indices. 



2 The Resistive Relativistic MHD Equations 

2.1 Conservative Formulation as a Stiff Hyperbolic Balance Law 

The total energy-momentum tensor of the system is made up by two contri- 
butions, T^ v = Tfif + Tj V . The first one is due to matter 

T^ = uuV+ PV ^, (2) 

where is the four velocity of the fluid, while u and p are the enthalpy and 
the pressure as measured in the co-moving frame of the fluid. The second 
contribution is due to the electromagnetic field 

Tf = F\F» X - {(F x «F XK )r)^ , (3) 

where F^ v is the electromagnetic tensor. If we introduce a laboratory observer 
defined by a four velocity = (1, 0, 0, 0), then the fluid four velocity w M and 
the standard three velocity in the laboratory frame are related as v = v % = 
u l /T, where T = (1 — v 2 ) -0 - 5 is the Lorentz factor of the fluid with respect to 
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the laboratory frame. The electromagnetic tensor, on the other hand, can be 
written as 

F v» = n ^E u -E^n u + e^ XK B x n K , (4) 

where E v and B v are the two spatial vectors (E° = B° = 0, E l = E i} B % = Bi) 
representing electric and magnetic field, respectively, while e M1/AK = [pp\K\ 
is the completely antisymmetric spacetime Levi-Civita tensor, with the con- 
vention that e 0123 = 1. The equations of motion can be derived from the 
conservation laws 

d,T^ = (5) 

and from the continuity equation 

dpipu") = 0, (6) 

where p is the rest mass density of the fluid. The electromagnetic field, on the 
other hand, obeys the Maxwell equations expressed in the form 

d^F*^ = 0, d v F^ = I", (7) 

where F*^ u = \e^ uXK F\ K is the dual of the electromagnetic tensor, while is 
the four vector of electric currents. The equations of resistive MHD differ from 
those of ideal MHD mainly because the second couple of Maxwell equations 
((TJ), accounting for the time evolution and for the divergence of the electric 
field, need to be explicitly solved. This requires that a relation is given between 
the currents and the electromagnetic field, the so called Ohm's law. In its 
most general form, the relativistic formulation of Ohm's law is a non linear 
propagation equation [20], but here, as in [22]. we will simply assume that 

I" = q u» + (tF^Uv, (8) 

where p c is the charge density in the co-moving frame while a is the electric 
conductivity. From ^ we easily derive the following expression for the spatial 
current vector 

J = Pc v + aT[E + vx B- (E -v)v\ , (9) 

where p c is the charge density in the laboratory frame. As done by [22J and 
[28] . to whom we address the reader for further details, we take care of the 
divergence-free character of the magnetic field by adopting the divergence 
cleaning approach presented in [7], namely by introducing two additional scalar 
fields \I/ and $ that propagate away the deviations of the divergences of the 
electric and of the magnetic fields from the values prescribed by Maxwell's 
equations. In total, the full set of RRMHD equations include the five equations 
for the fluid, plus the six equations for the evolution of the electric and of the 
magnetic field, plus the two equations about the divergences of the two fields, 
plus one more equation expressing the conservation of the total charge. In 
Cartesian coordinates, using the abbreviations d t = and = they can 
be written as: 
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d t D + diiDv*) = 0, 

d t Sj + d t z) = o, 

d t T + diS* = o, 

d t E l - e ijk djB k + d& = - J\ 

dtB i + e WdjE k + <9;$ = 0, 

d t * + diE' = Pc - 

<9 t p c + di J 1 = 0, 



where the conservative variables of the fluid are 
D = pT, 

S* = o;rV + e ijk E j B k , 
t = ujY 2 - P +UE 2 + B 2 ) , 



(10) 

(11) 

(12) 
(13) 
(14) 
(15) 
(16) 
(17) 

(18) 
(19) 
(20) 



expressing, respectively, the relativistic mass density, the momentum density 



and the total energy density. The spatial tensor Zj in (11), representing the 
momentum flux density, is 

" ~" " ' ' ' ' (21) 



Z) = ljT 2 v 1 Vj - E l Ej — B % Bj + p + \(E 2 + B 



where 5j is the Kronecker delta. In the rest of the paper we have assumed the 
equation of state of an ideal gas, namely 



p = (7- l)pe = 7i(w-p), 



(22) 



where 7 is the adiabatic index, 71 = (7— l)/7, e is the specific internal energy 
and uj = pe + p + p is the enthalpy. The system of equations ( 10 )-( 17) is 



written as a hyperbolic system of balance laws as in ([T]) , it has source terms in 
the three equations (13) that are potentially stiff (see [28J for a more detailed 
description of the different limits of the resistive MHD equations) and, as 
such, it can be treated with the procedure proposed by [TT], as we will show 
in Sect. [U 



2.2 Closed form recovering of the primitive variables from the conservative 
ones 



A fundamental difference with respect to the ideal MHD case is that the 
augmented set of conservative variables of the resistive relativistic equations 
allows for the recovering of the primitive variables from the conservative ones 
in closed form, at least when the equation of state is that of an ideal gas. This 
can be seen in the following way. Firstly, we shift the cross pruduct E x B 
from the right hand side to the left hand side of Eq. (19), then we square it, 
and we obtain 

(S-ExB) 2 = lu 2 T 2 (T 2 -1) . (23) 
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On the other hand, from (20) we obtain the enthalpy u as 



r 2 -71 



where we have used p = jx(co — p) as in (22). After replacing (24) into (23), 
simple calculations lead to the following quartic equation in the unknown 
Lorentz factor T as follows: 

A 4 r 4 + A 3 T 3 + A 2 T 2 + A ± T + A = 0, (25) 

where 

A i = C 1 -C*, A 3 = 2C 2ll D, A 2 = C 2 2 -2C lll - 1 2 D 2 , (26) 

A 1 = -2C 2ll D, A = 1 2 (C 1 + D 2 ), (27) 

with C x = (S - E X B) 2 , and C 2 = r - \{E 2 + B 2 ). The quartic S can 
be solved either analytically using the approach of Ferrari and Cardano [5] or 
numerically via a Newton-Raphson scheme. In our numerical experiments we 
found that for the purpose of accuracy and robustness, it is advisable to solve 
the quartic first analytically and then to improve the accuracy of the result 
by one or two additional Newton iterations. This is necessary since the com- 
putations of the roots for the analytical solution of the quartic may introduce 
a significant amount of roundoff errors on finite precision computer hardware 
even when using double precision arithmetic. It is only for this reason that the 
additional Newton iterations are performed. This step would not be necessary 
with exact arithmetic. 



As already pointed out by [36J, it turns out that (25) has two complex conju- 
gate solutions, plus two real solutions, of which only one is larger than unity, 
as physically required. Once the Lorentz factor is known, the other primitive 
variables can be computed in a straightforward manner. 



3 Numerical Method 



3.1 The PnPm Reconstruction Operator on Unstructured Meshes 



The main ingredient of the proposed numerical method to reach high order 
of accuracy in space is the PnPm reconstruction operator on unstructured 
meshes first introduced in [10J. It is a direct extension of the reconstruction 
algorithm proposed in ^2llT3] for finite volume schemes. For the details, we 
refer to the above mentioned publications and give only a short review in this 
section. The computational domain Q is discretized by conforming elements 
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Qi, indexed by a single mono-index i ranging from 1 to the total number of 
elements Ne- The elements are chosen to be triangles in 2D and tetrahedrons 
in 3D. The union of all elements is called the triangulation or tetrahedrization 
of the domain, respectively, 

N E 

Qn = {J Qi- (28) 
i=i 

At the beginning of a time-step, the numerical solution of ([I]) for the state 
vector W, denoted by Uh, is represented by piecewise polynomials of degree 
N from the space 14, spanned by the basis functions = $i(x), i.e. at t = t n 
we have for each element 

UfcfoO = 5>,(£K. (29) 
I 

From the polynomials Uh, we then reconstruct piecewise polynomials Wh of 
degree M > N from the space Wh, spanned by the basis functions = ^i(x): 

w h (x,t n ) = Y,Vi(x)ti?. (30) 
i 

As stated in [10], the form an orthogonal basis and are identical with the 
up to polynomial degree N. For performing the reconstruction on element 
Qi, we now choose a reconstruction stencil 

Si = U Qm (si) 

k=l 

that contains a total number of n e elements. Here 1 < k < n e is a local 
index, counting the elements in the stencil, and j = j(k) is the mapping from 
the local index k to the global indexation of the elements in Qq. For ease of 
notation, we write in the following only j, meaning j = j(k). 

In the present paper the three operators 

t n+l 

(/, g) Ql = J J (m t) ■ g(x, *)) dv dt, (32) 

[/,^ = / (f(x,t).g(x,t))dV, (33) 

f n + l 

{/.»W= / / (f(x,t)-g(x,t))dSdt, (34) 

denote the scalar products of two functions / and g over the space-time el- 
ement Qi x [t n ;t n+1 ], over the spatial element Qi, and over the space-time 
boundary element dQi x [t n ; t n+1 ] respectively. The operators (/, g) and [/, gf, 
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written without the index Q i} will denote scalar products on the space-time 
reference element Qe x [0; 1] and on the spatial reference element Qe at time 
t, respectively. The spatial reference element Qe is defined as the unit simplex 
with vertices (0, 0), (1, 0), (0, 1) in two space dimensions and vertices (0, 0, 0), 
(1, 0, 0), (0, 1, 0) and (0, 0, 1) in three space dimensions, respectively. 

The reconstruction is now obtained via /^-projection of the (unknown) piece- 
wise polynomials Wh from the space Wh into the space Vh on each stencil Si 
as follows: 

[^Mq^I^Mq^ VQjeS t . (35) 



Note that during the reconstruction step, the polynomials Wh are continuously 
extended over the whole stencil <Sj. After reconstruction, the piecewise poly- 
nomials Wh are again restricted onto each element Qi. The number of elements 
in the stencils are chosen in such a way that the number of equations in (35) 
is larger that the number of degrees of freedom in the space Wh- Eqn. (35) 
constitutes thus an overdetermined linear algebraic equation system for the 
coefficients of Wh and is solved using a constrained least squares technique, see 
[TUfT2] . The linear constraint is that Eqn. (35) is at least exactly satisfied for 
Qj = Qi, i.e. inside the element Qi under consideration. The integral on the 
left hand side in (35) is computed using classical multidimensional Gaussian 



quadrature of appropriate order, see [31] . The integral on the right hand side 
can be computed analytically and involves the standard element mass-matrix. 
The resulting M-exact PnPm least squares reconstruction can be interpreted 
as a generalization of the A;-exact reconstruction proposed for pure finite vol- 
ume schemes by Barth and Frederickson in their pioneering work [I]. 



3.2 The Local Space- Time Discontinuous Galerkin Predictor for Stiff Balance 
Laws 



The original ENO scheme of Harten et al. [IT] as well as the ADER-FV and 
ADER-DG schemes developed by Titarev and Toro [32] and Dumbser and 
Munz [H] use the governing PDE itself in its strong differential form to ob- 
tain high order of accuracy in time. This is achieved via the so-called Cauchy- 
Kovalewski procedure that substitutes time derivatives with space derivatives 
via successive differentiation of the governing PDE with respect to space and 
time. This procedure becomes very cumbersome or even impossible for gen- 
eral nonlinear hyperbolic PDE systems. In [TT], [TO] and [3] an fully numerical 
approach was presented that replaces the semi-analytical Cauchy-Kovalewski 
procedure by a local weak formulation of the governing PDE in space-time. 
While the approach presented in [11] relies on a local discontinuous Galerkin 
approach in space-time, which is also able to handle stiff source terms, the 
methods proposed in [TD] uses a local continuous Galerkin method in space- 
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time as predictor. In the present article we use the local space-time discontin- 
uous Galerkin method due to the stiffness of the source terms. 
We underline that the local space-time DG method is only used as local pre- 
dictor for the construction of a solution of the PDE in the small, as it was 
called by Harten et al. in [17J. The local space-time predictors are then in- 
serted into a global corrector, which is fully explicit and provides the coupling 
between neighbor cells. As a consequence, the resulting nonlinear algebraic 
systems of the local space-time Galerkin methods are element local and not 
globally coupled, as in the global and fully implicit space-time Galerkin ap- 
proach introduced by van der Vegt and van der Ven [3"4"1T21~] . 



The details of the local space-time DG predictor method are already given 
in [UJ, hence we will only briefly recall the basic ideas here. We start from 
the strong formulation of PDE Q and transform the PDE into the reference 
coordinate system (£, r) of the space-time reference element Qe x [0; 1] with 
£ = (£, rf) and being the nabla operator in the £ — rj reference system: 

d 

— W + V c • F* (W) = S*. (36) 
The modified flux tensor and the modified source term are given by 

F* := AtF(W)J T , S* := AtS(W), J = (37) 

ox 



as revealed by simple algebraic manipulations. We now multiply Eqn. (36 ) by a 
space-time test function 9 k = rj, r) from the space of piecewise space-time 
polynomials of degree M and integrate over the space-time reference control 
volume Qe x [0; 1] to obtain the following weak formulation: 

U k , JUvA + (6 k , V € ■ F* (W h )> = (6 k , S* h (W h )) . (38) 



In the numerical solution of Eqn. (38) we use the same ansatz for as well 



as for the flux tensor and the source term , i.e. 

Wh = v, r) = £ 0z(£, V, := em, (39) 

i 

F* h = F*(£, V ,t) = Y: rj, r)F* := d{F*, (40) 
l 

S* h = S* h (Z, V ,T) = J2m,V,T)St -.= 6$ ■ (41) 
i 

The degrees of freedom of the flux F^ and the source S£ can be computed 
from the ones of the state vector Wj either via the more accurate but also 
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more expensive L 2 -projection, 



,0,)F? 



«,F*(W ft )), (Mi)Sf = (0*,5*(W fc )>, 



(42) 



or in a simple and cheap nodal fashion, if a nodal space-time basis as the one 
in [TO] is used: 



F*(Wi), 5,* = 5*(Wj; 



(43) 



In the element-local weak formulation of Eq. (38) we apply integration by 
parts to the first term and thus obtain 



KWh} 1 -[9 h ,w h ] 



d_ 
dr 



9 k ,W h ) + (9 k ,V r F* h ) = (9 k ,S* h ), 



(44) 



where the initial condition at relative time r 
weak sense by the term [9 k , Whf 



is taken into account in a 
We recall that Wh is the piecewise polynomial 



obtained by the high order P^Pm reconstruction operator summarized in |3.1 



We also note that the first two terms in (44) correspond to the choice of 



an upwind flux in the time direction, which is consistent with the causality 
principle that states that no effect can occur before its cause. Inserting the 



ansatz (39)-(41) into (44) we obtain 



' d 

dr 



■ Oi) SI . (45) 



After defining the following universal matrices (that need to be computed only 
once on the reference element) K x = [9 k ,9\\ l — (^J^, 6^, K 

M=(6 k ,6 l ),F 
notation: 



(9 k ,'V^9i), 

[9 k ,^i}° we can rewrite (45) in the more compact matrix 



KxW, + K r F* t = F w? + M5;. 



(46) 



Eqn. (46) is an element-local nonlinear algebraic system for the unknowns Wj. 



For its solution we use the following simple iterative scheme, similar to the 
one proposed in [TO] : 



(K 1 )- 1 F u ? r - (Ki)" 1 ^ • FJ 



(47) 



As in [10] the matrices contained in (K x ) have the remarkable prop- 



erty that all their eigenvalues are zero, which makes (47) a contractive fixed 
point iteration in the homogeneous case (i.e. when S — 0) and thus existence, 
uniqueness and convergence to the unique solution are guaranteed by the Ba- 
nach fixed point theorem. Furthermore, in the linear case, the method is even 
guaranteed to converge to the exact solution in M + 1 steps from any initial 
guess. In the non-homogeneous case with stiff source terms, however, it is nec- 
essary to take the source implicitly, which is done in the present paper. We 
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use the following simplified model for the implicit source term: 



<?*,i . . 9S f^i+i tT7 4\ dS dS dV dV ( dW 



Sl ' * s i' +At dw\ Wl ~ Wl J> dW dVdW' dW \dV 

'(48) 

where the Jacobian of the source with respect to the conservative variables W 
is computed by the chain rule, taking first the derivatives with respect to the 
vector of primitive variables V. The derivative of V with respect to W can 
be computed easily by the theorem on the derivative of the inverse function. 
To simplify the computations, we evaluate only once per iteration at the 
current space-time average value of Wh- 
in our numerical experiments we also found that in the very stiff case, the 
choice of the initial guess W; seems to be very crucial. We therefore adopt the 



following strategy: First we solve (47) at the first order level, which becomes 



a simple Newton-Raphson scheme for the space-time cell-average W as 

/ = W - S* (W) - v% = 0, (49) 



where the initial guess W° = u" is used for all variables apart from the electric 
field. For the electric field we use E obtained from the relaxation of E and 
v to equilibrium assuming the stiff limit o — > oo and holding all the other 
conservative variables in W constant. In our experiments the Newton method 



applied to eqn. (49) with this initial guess typically converges to machine zero 
(10 -14 ) after two or three iterations and is robust even for very large values of 
a, such as a = 10 12 . The resulting cell average W is then used as initial guess 



for the high order space-time solution of eqn. (47), i.e. we set Wjj = W. 



3.3 The Fully Discrete PnPm Schemes 



The fully discrete one-step form of the proposed PnPm schemes is derived as 
follows: we first apply the operator -)q, to PDE Q and obtain 

+ V • F(W)) — S(W)) . (50) 




For the first term in Eqn. (50) we approximate W with Uh from the space Vh 
and perform integration by parts in time. Note that the do not depend on 
time and therefore their time derivatives vanish. For all the other terms in Eqn. 



(50) the vector W is approximated by the solution Wh of the local space-time 



discontinuous Galerkin predictor of section |3.2| Since W4 will usually exhibit 
jumps at the element boundaries, we introduce a numerical flux to resolve 
these jumps. We hence obtain the following family of fully discrete one-step 
P N P M schemes for PDE <Q: 
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,n + l 



+ {^ k ,G i+h (W^,W^-n} aQi = (* k ,S(yV h )) Qi 



(51) 



where denotes the boundary extrapolated data from within element Qi 
and denotes the boundary extrapolated data from the neighbor, respec- 
tively. In the test sections of this paper, we use the Rusanov flux for 
which in the case of the resistive relativistic MHD equations becomes par- 
ticularly simple. Due to the presence of the full Maxwell equations, whose 
maximum eigenvalue is the speed of light, i.e. A max = 1, it reduces to 



G l+h (W„ W+) -n=\ (F(W+) + F(W h 



1 



n 



(52) 



As an alternative, we also propose the following strategy, which gives slightly 
better results for the hydrodynamic quantities: For the evolution of the hy- 



drodynamics, Eq. (10 )-( 12 ), one can use the HLL flux 



G 4+| (W,,W+)-n 
with 



(a+F(W+) + a-F(W ft " 



• n — a + a" 



<v + a 



,+ 



max{0, \f , \~f }, 



max{0, -X s , -A+}, 



(53) 
(54) 



where Xf and X s denote the fastest and the slowest of the ideal MHD mag- 
netosonic speeds along the direction of the flux, and computed through the 
exact or approximate solution of the corresponding quartic as in [8] . 

For a quadrature-free implementation that requires only the solution of one 
Riemann problem per space-time element interface we refer the reader to [10] • 



4 Numerical Test Cases 



In this section we present some of the test cases of Palenzuela et al. [28J who 
used a second order accurate TVD scheme with IMEX Runge-Kutta time- 
integration on Cartesian meshes. In the rest of the section we use schemes of 
order better than two in space and time on unstructured triangular meshes 
and the constant k in the equations (15) and (16) for the divergence cleaning 
is set equal to unity in all tests. 
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4-1 Large Amplitude Alfven Wave 



This smooth unsteady test case with exact analytical solution was introduced 
for the ideal relativistic MHD equations by Del Zanna et al. jH] and was solved 
for the first time on unstructured triangular meshes with high order P^Pm 
schemes in pU]. Since the resistive MHD equations tend asymptotically to the 
ideal ones in the stiff limit (a — > oo), this is an ideal test case to assess the 
accuracy of our scheme in the stiff limit of the governing PDE system. 
The test case consists of a periodic Alfven wave whose initial condition at 
t = is chosen to be p = p = 1, B = B (1, cos (kx) , sin (kx)) T , v = —va/B ■ 
(0, B y , B Z ) T , E = —v x B and <fi = ip = q = 0. We furthermore use the 
parameters k = 2n, 7 = § and B = 1, hence the advection speed of the 
Alfven wave in x-direction is va = 0.38196601125 (see [S] for a closed analytical 
expression for va)- The 2D computational domain is Q = [0; 1] x [0; 0.4] with 
four periodic boundary conditions, and the final time corresponding to an 
entire advection period is t — 1/va = 2.618033988. The initial condition 
represents the exact reference solution to be compared with our numerical one. 
Since this test case was constructed for the ideal relativistic MHD equations, 
we have to use a rather stiff value for the conductivity (a = 10 7 ) in the resistive 
case to reproduce the ideal equations asymptotically. For the fifth order P\P± 
scheme this has shown to be not enough to get the full order of accuracy, 
hence in this case we even use a = 10 8 . In all our computations a constant 
Courant number of CFL = 0.5/ (2N + 1) is used. 

A representative unstructured triangular mesh is visible in the left panel of 
Fig. [T] together with a surface plot of the quantity B z . In the right panel we 
compare the exact solution after one period with the numerical one obtained 
on a very coarse mesh of 8 triangles on the x-axis using the -P1-P4 scheme. For 
this purpose, the reconstructed fourth degree polynomials are evaluated at the 
final time on 100 equidistant points along the x-axis in order to make use of 
the high order polynomial sub-cell resolution contained in each element. We 
emphasize the excellent agreement with the exact solution even on this very 
coarse mesh. Note that with the TVD scheme used in [28] there were clearly 
visible errors even on a fine mesh using 50 points along the rc-axis. 
Table [l] shows the errors and the orders of convergence measured in the L 2 
norm for the flow variable B y . The number Nq denotes the number of triangle 
edges along the x-axis. We stress that the P1P4 scheme on the very coarse 
mesh with Nq = 8 allows to achieve an accuracy higher than the P0P2 scheme 
on the finest mesh with Nq = 64. The nominal order of accuracy M + 1 
has been reached for all PnPm schemes under consideration. In it was 
reported that when using IMEX Runge-Kutta schemes for time-discretization 
the authors encountered problems with the convergence rates for the relaxed 
variables, i.e. for the electric field that suffers from the presence of the stiff 
source terms. With our local space-time Galerkin predictor method, where 
nonlinear flux and source term are fully coupled in the predictor stage and 
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Table 1 

Large Amplitude Alfven Wave. Convergence study of P^Pm schemes from third to 
fifth order of accuracy, a = 10 7 , apart from the P\Pt± scheme where a = 10 8 . Errors 
are computed for variable B y . 



P0P2 


P1P2 


P2P2 


N G 


L 2 




N G 


L 2 




N G 


L 2 




16 


1.71E-02 




8 


9.12E-04 




8 


8.97E-04 




24 


5.32E-03 


2.9 


12 


2.26E-04 


3.4 


12 


2.92E-04 


2.8 


32 


2.26E-03 


3.0 


16 


9.34E-05 


3.1 


16 


1.67E-04 


1.9 


64 


2.79E-04 


3.0 


24 


2.53E-05 


3.2 


24 


4.98E-05 


3.0 


P0P3 


P1P3 


PiPa 


N G 


L 2 




N G 


L 2 


L 2 


N G 


L 2 




12 


1.81E-03 




4 


7.18E-03 




4 


3.32E-03 




16 


4.52E-04 


4.8 


8 


3.75E-04 


4.3 


8 


2.95E-05 


6.8 


24 


7.35E-05 


4.5 


12 


7.91E-05 


3.8 


12 


4.46E-06 


4.7 


32 


1.98E-05 


4.6 


16 


2.82E-05 


3.6 


16 


1.07E-06 


5.0 



where the optimal local space-time polynomial distribution is found due to the 
Galerkin orthogonality property, such problems have not been encountered. 
Therefore and for the sake of completeness, we show the convergence rates 
for the relaxed variable E y in Table [2] for the schemes P0P2, P0P3 and P1P4. 
We deduce from the results of Table [2] that the nominal order of accuracy is 
reached even for the electric field, which contains the stiff source term. This 
confirms the results already presented in [TT] , where uniform order of accuracy 
in space and time was found in the stiff as well as in the non-stiff case. 




Fig. 1. Large Amplitude Alfven Wave. Left panel: very coarse unstructured trian- 
gular mesh (h = 1/8) used for the fifth order P1P4 scheme and surface plot of the 
quantity B z at the final time t = 2.618033988. Right panel: Comparison of exact and 
numerical solutions for B y at the final time obtained with the P\Pa scheme on the 
very coarse mesh. A cut along the line y = is shown, evaluating the reconstructed 
polynomials on 100 equidistant points. 
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Table 2 

Large Amplitude Alfven Wave. Verification of the order of accuracy for a variable 
affected by the stiff source term. We use the quantity E y and some selected PnPm 
schemes. 



P0P2 


PqPs 


PiPl 


N G 


L 2 




N G 


L 2 


<D L 2 


N G 


L 2 




16 


7.66E-03 




12 


6.09E-04 




4 


9.43E-04 




24 


1.90E-03 


3.4 


16 


2.11E-04 


3.7 


8 


1.22E-05 


6.3 


32 


7.75E-04 


3.1 


24 


4.02E-05 


4.1 


12 


2.06E-06 


4.4 


64 


9.56E-05 


3.0 


32 


1.14E-05 


4.4 


16 


5.18E-07 


4.8 



4-2 Self-similar Current Sheet 



This smooth test case was first proposed by Komissarov et al. [23J , it has been 
presented also in Palenzuela et al. [28J and it provides a truely resistive test, 
far from the ideal MHD limit. It has the following exact analytical solution 
for the y-component of the magnetic field: 

B y (x,t) = B erf(±^x y ) , (55) 

where erf is the error function. The initial time for this test case is t — 1 
and the initial condition at t — 1 is given by p — 1, p — 50, E — v — 
and B = (0, B y (x, 1), 0) T . We choose 7 = | and B = 1. The conductivity 
is chosen as a = 100, which means a moderate resistivity. The problem is 
solved with two different fourth order PnPm schemes on the two-dimensional 
computational domain Q = [—1.5; 1.5] x [—0.5; 0.5], where we impose periodic 
boundary conditions in ^-direction and Dirichlet boundary conditions consis- 
tent with the initial condition in x-direction. The first scheme is a pure finite 
volume method (P Q P 3 ) using the component- wise WENO reconstruction pro- 
posed in [12], running on a mesh with h = 3/32, which corresponds to an 
equivalent one-dimensional resolution of 32 points. The second scheme is the 
P2-P3 method which is part of the new intermediate class of numerical schemes 
discovered in [10], running on a very coarse mesh with h = 3/8, i.e. using only 
8 points in the one-dimensional case. The mesh is depicted on the left panel of 
Fig. [2] together with a surface plot of the magnetic field in ^/-direction. Both 
numerical solutions are compared at time t = 10 with the exact solution given 
by Komissarov et al. [23] and Palenzuela et al. [28] on the right panel of Fig. 
[2] We note an excellent agreement with the exact solution and underline that 
the use of high order methods in space and time allows us to use very coarse 
meshes, compared to standard second order TVD schemes. 
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Fig. 2. Self-similar current sheet. Left panel: unstructured triangular mesh used for 
the -P2-P3 scheme and surface plot of the quantity B y . Right panel: Comparison of 
exact and numerical solutions at time t = 10 obtained with two different fourth 
order P^Pm schemes on different meshes. A cut along the line y = is shown, 
evaluating the reconstructed polynomials on 100 equidistant points. 

Table 3 

Initial states left (L) and right (R) for the relativistic MHD shock tube problem. 
The last column reports the final time t e considered in the numerical test. 



Case 


P 


P 


u 


V 


w 


By 


B z 


B x t e 


L 


1.08 


0.95 


0.4 


0.3 


0.2 


0.3 


0.3 


2.0 0.55 


R 


1.0 


1.0 


-0.45 


-0.2 


0.2 


-0.7 


0.5 


2.0 



4-3 Shock Tube Problems 



In this section we solve the fifth of a series of test problems proposed by 
Balsara in [TJ. We solve the RRMHD equations with different values for the 
conductivity a. The initial condition is given by two piecewise constant states 
separated by a discontinuity at x = 0. The left and right values for the prim- 
itive variables are reported in Table [3 Furthermore, we set E = —v x B, 
cj) = ip = q = and 7 = |. The conductivities in our test cases are chosen as 
a — 0, a — 1, a — 10, a — 10 2 , o = 10 3 and a = 10 6 . The computational do- 
main is Q = [—0.5; 0.5] x [0; 0.05] with periodic boundaries in y-direction and 
Dirichlet boundaries consistent with the initial condition in x-direction. We 
use an unstructured triangular mesh of characteristic size h = 1/400, which is 
depicted together with a surface plot of the density p in Fig. |3j A cut through 
the solution along the x-axis is shown in Fig. [4] for all different values of a 
used in this series of test cases. The exact solution is the one for the ideal 
RMHD equations, as published in [15]. The essential wave structures of the 
ideal RMHD Riemann problem can be noted for a = 10 3 or greater. For values 
below, the resistivity leads to a considerable diffusion of the discontinuities. 
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Fig. 3. RRMHD shock tube test problem using a P0P2 WENO scheme and different 
values for the conductivity a. The unstructured triangular mesh is shown together 
with a surface plot of the density p. Left panel: a = 10 3 . Middle panel: a = 10 2 . 
Right panel: a = 10. A cut along the line y = is shown at 400 equidistant points. 



4-4 Rotor Problem 



In this section we solve a resistive relativistic version of the MHD rotor prob- 
lem proposed by Balsara and Spicer [2] . Our computational setup is a variation 
of the ideal relativistic MHD rotor test case of Del Zanna et al. [35J. In con- 
trast to [35J, who solved the ideal RMHD equations on a perfectly regular 
Cartesian mesh, we solve this test case in Cartesian coordinates on a circular 
computational domain with radius R — 0.5 using an unstructured triangu- 
lar mesh with a characteristic mesh spacing of h = 0.004 towards the center 
and h = 0.005 at the outer border of the domain, leading to a total number 
of 72320 triangles. The rotor has a radius of i?o — 0.1 and is spinning with 
an angular frequency of io s = 8.5, leading to a maximal toroidal velocity of 
= 0.85. The density is p = 10 inside the rotor and p = 1 in the outer fluid 
at rest. The pressure is p — 1 and the magnetic field is B — (1, 0, 0) T in the 
whole domain. The initial electric field is, as usual, computed as E = —v x B. 
We use a P0P2 scheme with component-wise WENO reconstruction. No taper 
is applied to the initial condition, as in [35J, and 7 = 4/3 is used. Transmissive 
boundary conditions are applied at the outer boundaries. The computational 
domain and the results for the pressure at time t = 0.3 are shown in Fig. |5]for 
different values of the electric conductivity. We solve the problem with a = 10 
and a = 10 5 and, as a reference solution, we also show the results obtained 
with the ideal RMHD equations. The ideal RMHD results agree qualitatively 
very well with those obtained with the RRMHD equations using the larger 
conductivity a = 10 5 . For the case of a lower conductivity (a = 10) one can 
clearly see that the wave structure is completely different, with a faster mov- 
ing electric field that is governed directly by the Maxwell equations and no 
longer resulting from the relation E = —v x B as in the ideal case. 
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Fig. 4. RRMHD shock tube test problem using a P0-P2 WENO scheme and different 
values for the conductivity cr. The exact solution is shown for the ideal RMHD 
equations. The density p is plotted on the top of the figure and the magnetic field 
component B y on the bottom. 
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4-5 Orszag-Tang Vortex 



The last of this series of test cases is a resistive relativistic analogous of the 
Orszag-Tang vortex problem [27] studied extensively in [6|29II19] . The compu- 
tational domain is Q = [0; 27r] 2 . The initial condition of the problem is given 
by 

(p,u,v,p,B x ,B y ) = (1, -sin(y),sin(x), 1, -B sin(y), B sin(2x)) , (56) 

with w = B z = and 7 = |. The problem is solved up to t = 4.5 using a 
P0-P2 scheme with component-wise WENO reconstruction on an unstructured 
triangular mesh with 55292 elements (h = ^). The results are shown for the 
pressure in Fig. [6] for times t = 0.5, t = 2.0 and t = 4.5 using two different 
values for the conductivity, a — 10 and a = 10 3 . As in the original problem 
[27] , the smooth sinusoidal initial condition evolves in time to form complex 
shock dominated structures for the large value of the conductivity. For small 
conductivities, much less waves are present due to the diffusion caused by the 
electric resistivity. 



5 Conclusions 



In this paper we have solved the resistive relativistic magnetohydrodynamics 
equations using the class of methods introduced in Dumbser et al. JT0j and 
named PnPm schemes. The equations present source terms that are poten- 
tially stiff when the ideal limit of infinite conductivity is recovered. As such, 
they are naturally accounted for through the application of the local space- 
time discontinuous Galerkin predictor, originally deviced in [11]. 

To our knowledge, the computations presented in this paper are the first better 
than second order accurate simulations in space and time ever done for the stiff 
limit of the RRMHD equations and the results obtained point to favour higher 
order methods over standard second order TVD schemes. In particular, the 
accuracy that can be achieved with high order P^Pm schemes on very coarse 
meshes makes them promising tools for simulations of physical processes that 
require high computational resources, such as a large class of time dependent 
problems involving magnetic reconnection in astrophysical context. Further 
directions of future improvement are represented by the generalization of the 
scheme into full general relativity as well as the inclusion of more complex 
Ohm's laws. 
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Fig. 6. Pressure field for the resistive relativistic Orszag-Tang vortex problem at 
times t = 0.5, t = 2.0 and t = 4.5 (from top to bottom). Left column: a = 10. Right 
column: a = 10 3 . 
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